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ABSTRACT 


In  an  effort  to  understand  certain  ideas  and  concepts  associated  with 

fft  C  ‘ 

multi-grid  iterations  we"  give  an  in-depth  study  of  a  particular  simple 

l  *  ' 

problem.  We'consider  a  standard  finite-difference  system  associated  with 
—the-  two-point  boundary  value  problem t 


-(pu1)’  +  bu1  +  qu  =  0,  u(0)  =  u(l)  =  0  . 


9  h  h 

The  operators  1^  ,  l£h  are  "operator"  based  interpolation  and  pro¬ 
jection  operators  while  the  smoothers  are  the  damped  Jacobi  iterations  with 
parameter  a  >  0  . 

We  determine^the  exact  rates  of  convergence  for  the  "two-grid"  scheme 
and  upper  bounds  (  <1  !  )  for  the  multi-grid  schemes.  Experimental  results 


are  discussed. 


1.  Introduction 


The  multi-grid  approach  for  the  numerical  solution  of  boundary  value  prob 
lems  for  elliptic  partial  differential  equation  is  proving  itself  as  one  of 
the  fastest  and  most  efficient  methods  -  see  [1],  [3],  [4],  [5],  [6],  [17]. 
Moreover,  there  are  a  large  number  of  theoretical  papers  on  this  subject  - 
see  [2],  [3],  [6],  [7],  [8],  [9],  [10],  [11],  [13],  [16].  Nevertheless,  it 
seems  (at  least  it  seems  so  to  these  authors)  that  we  are  just  beginning  to 
understand  this  powerful  idea.  In  particular,  there  are  questions  of:  how 
do  we  choose  the  interpolation  and  projection  operators?,  how  do  we  choose 
the  smoothing  operators?,  what  do  we  mean  when  we  say  smoothing?  and  ...? 

This  report  is  a  reflection  of  our  efforts  to  understand  and  appreciate 
the  theoretical  insights  of  Frederickson  [7],  McCormick  and  Ruge  [13], 
McCormick  [14],  [15]  and  Greenbaum  [8]  and  apply  those  ideas  to  extend  the 
explicit  convergence  rates  given  by  Hackbusch  [12,  (2.21)],  [11]  for  the  very 
simplest  problem 


u"  =  f,  u(0)  =  u(l)  =  0  . 

Specifically,  we  consider  the  two-point  boundary-value  problem 

(1.1)  Lu:  =  -(pu')'  +  b(x)u'  +  qu  =  f,  0  <  x  <  1 

(1.2)  u  ( 0 )  =  u(l)  =  0 
where  p(x),  b(x),  q(x)  are  smooth  functions  and 

(1.3)  p(x)  >  P0  >  0,  q(x)  >  0  . 


In  section  2  we  describe  a  basic  approach  to  multi-grid  which  is  based  on 
the  ideas  of  Frederickson  [7],  McCormick  and  Ruge  [13]  and  Greenbaum  [8]. 
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In  section  3  we  describe  a  discretization  (finite-difference)  of  the 

problem  (1.1),  (1.2)  and  a  specific  two-grid  iterative  procedure  for 

its  solution.  In  section  4  we  describe  the  class  of  damped1  Jacobi 

"smoothers"  and  use  our  knowledge  of  these  schemes  and  their  eigenvectors 

h  2h 

to  describe  the  basic  spaces:  Range  l£h  =  R  and  Nullspace  1^  Lh  =  n  • 

In  section  5  we  obtain  estimates  for  the  norm  decay  of  a  single  step  in  a 
two  grid  scheme  for  two  different  norms.  In  addition  we  obtain  a  better 
estimate  for  the  norm  decay  for  all  iterative  steps  beyond  the  first  in 
both  these  norms.  Interestingly  enough,  this  estimate  is  the  same  in  both 
norms.  This  latter  result  is  an  improvement  over  the  estimates  of  Hackbusch 
[11],  [12]. 

It  is  well-known  that  the  problem  (1.1),  (1.2)  is  equivalent  to  a 
self-adjoint  problem.  Moreover,  the  discretization  (3.4)  is  also  equivalent 
to  a  symmetric  problem.  In  fact,  our  multi-grid  treatment  of  this  problem 
is  equivalent  to  the  "same"  multi-grid  treatment  of  this  symmetric  problem. 
This  equivalence  is  not  needed  for  the  discussion  in  sections  1-5.  However, 
as  we  turn  to  the  extension  of  a  two-grid  scheme  to  a  true  multi-grid  scheme, 
we  require  this  information.  In  section  6  we  demonstrate  this  equivalence. 

In  section  7  we  describe  the  n-grid  "saw-tooth"  multi-grid  schemes  and  give 
a  theory  (closely  related  to  a  theorem  of  McCormick  [14],  [15])  which  de¬ 
scribes  the  rates  of  convergence  of  this  scheme.  In  addition  section  7  con¬ 
tains  some  experimental  results. 


2.  A  Basic  Theory 

The  theory  presented  in  this  section  is  based  on  the  work  of 
Frederickson  [7],  McCormick  and  Ruge  [13]  and  Greenbaum  [8],  We  consider 
a  finite-dimensional  linear  space  and  a  problem 


(2.1) 


where 


L^U  =  f;  u,  f  e  Sh 


4r  sh  *  sh 


is  a  linear,  nonsingular  operator. 

Multi-grid  is  an  iterative  method  for  the  solution  of  this  problem. 
The  basic  idea  is  to  utilize  another  finite  dimensional  space  with 


(2.2) 


dim  S2h  <  dim  . 


Ph  u 

Hence  we  require  operators  ijj  ,  I^h  which  enable  us  to  effect  com¬ 
munication  between  these  spaces.  In  particular,  we  have 


(2.3a) 


(2.3b) 


I?*.  5  +  s 

xh  *  ^h  52h 


I  •  S  -*■  S 
x2h*  *2h  3h 


(projection) 


(Interpolation) 


2  h  h 

where  1^  ,  I^  are  linear  operators.  We  also  require  a  "smoothing" 
operator  and  a  "coarse  grid"  operator  L,,^.  The  smoothing  operator 
is  an  affine  operator  which  has  U,  the  unique  solution  of  (2.1)  as 
its  only  fixed  point.  That  is 


a 


where  G^:  -*•  is  a  linear  operator  and  if  U  is  the  solution  of 

(2.1),  then 


(2.4b) 


shu  =  u 


Finally,  the  "coarse  grid"  operator 


(2.5) 


L  •  S  S 

L2h‘  i2h  ^2h 


is  a  linear,  nonsingular  operator  taking  S2h  onto  itself. 
Let  U°  €  be  a  guess  for  the  solution  U  of  (2.1). 


(2.6a) 


c°  -  U  -  U°  , 


(2.6b) 


U  -  V”  , 


(2.6c) 


e  «  U  -  0  =  Gh(U-U°)  =  Ghc°  , 


(2.6d) 


r  =  f  -  LhU  =  Lh(U-U)  =  Lhe  , 


(2.6e) 


R  *  Ihv  • 


Solve 


(2.7) 


L2h£  =  R  »  1*e*»  e  =  L2hR 


(2.8) 


U  +  ,  e1  =  e 


11°:  =  U1 


and  return  to  (2.6a). 


r 


Remark:  While  one  might  say  that  we  have  merely  described  a  two-grid 
scheme,  the  iterative  scheme  described  above  does,  in  fact,  describe 
"multi-grid"  schemes.  The  point  is  that  the  "coarse  grid  correction" 
operator  may  be  a  complicated  procedure  involving  more  grids. 

The  work  of  Frederickson  [7],  McCormick  and  Ruge  [13]  and  Greenbaum 
[8]  suggests  we  study  two  fundamental  subspaces.  These  are 


(2.9a) 

(2.9b) 


R  =  range  of  l!jh  , 

?  h 

n  =  nullspace  of  . 


In  addition  we  consider  a  special  two  grid  "coarse  grid"  operator 

(2J°)  L2h  =  ljj*'LhI2h  • 


This  particular  operator  is  the  "Galerkin  choice"  and  is  "optimal"  in  a 
certain  sense.  This  fact  is  emphasized  by  the  following 

Lemma  2.1 :  Consider  one  iteration  as  described  above  by  (2.6a)-(2.8)  with 


(2.11) 


Suppose  e  e  R  .  That  is,  suppose  there  is  a  w(2h)  e  and 
(2.12)  e  =  e(h)  =  l!jhw(2h)  . 


Then 


(2.13a) 


i  =  w(2h) 


Hence,  the  problem  is  solved. 


Proof:  From  (2.6d)  and  (2.12)  we  have 

r  =  LhS  LhI2hw^2h^  * 

Thus 

R  ■  ■h'v  *  <ih\i5h>»<2h'  • 

That  is, 

(2.14)  R  =  L2hw(2h)  . 

Hence,  from  (2.7)  we  have  (2.13a).  Finally,  (2.13b)  follows  from  (2.8) 
and  (2.12).  | 

Now  suppose  we  can  write  Sh  as  the  direct  sum  (not  necessarily  orthog 
h  2h 

onal )  of  Range  (I2^)  anc*  Nulls-pace  (I^nL^).  That  is,  every  grid  function 
w(h)  e  Sh  can  be  uniquely  written  as 

(2.15a)  w(h)  =  I2hw](2h)  +  w2(h) 

where 

(2.15b)  I2hLhw2(h)  =  0  . 

Let  us  apply  this  decomposition  to  the  intermediate  error  e  =  e(h)  .  Then 


(2.16) 


?(h)  =  l"hw‘ (2h)  +  e'(h)  . 


Thus 

r  =  Lhl5hw](2h)  +  L^1  (h) 
and 


R  -  >h'V  *  • 

Using  lemma  2.1  we  see  that 

(2.17)  U1  =  0  +  l2hw1(2h) 


and  hence 


e1  =  U  -  U1  =  U  -  0  -  l2hw](2h)  . 


That  is 


1 

€ 


e1 (h)  . 


Thus,  the  convergence  of  the  process  can  be  tested  by 

ll«1(h)ll/l|c°ll 

where 

e^(h)  e  Nullspaae(l^\.^)  =  n 

and  the  norm  involved  is  any  norm.  The  authors  mentioned  above  all  dealt 
with  self-adjoint  Lh  and  assume  that 


In  this  case  it  is  an  easy  matter  to  see  that  the  Range  (I,.  )  and 

d  n 

9  h 

Nullspaoe  (I^nLh)  are  L^-orthogonal  complements  of  Sh  .  In  the  general 
case  we  must  assume  the  decomposition  (2.15a).  However,  it  is  a  simple 
counting  argument  to  see  that  (2.15a)  is  valid  if  L^h  is  non-singular  and 
i)  rank  l£h  =  dim  S2h 

ph 

ii)  dim  Nullspaoe  1^  dim  -  dim  . 

Recall  that  Lemma  2.1  implies  that  the  zero  vector  is  the  only  vector  common 
to  both  subspaces. 
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o 


a 


3.  The  Discrete  Problem 


Let  an  integer  N  >  1  be  chosen  and  set 


(3.1) 


2(N+1)  =  (2N+1 )+l 


and  let  the  fine  "grid  points"  be  given  by 


(3.2) 


x  . ( h )  =  jh,  j  =  0,1 ,2,...,2N+2  . 

J 


We  define  a  difference  operator  by 


(3.3a) 


where 


(3.3b) 


[Lh^k  ■  -Vk-1  *  Vk  -  >kuk+i 


’ 

„  -  Pk-^  +  bk 

ak  -  “Jr +  2h  * 


P*-ii  ,  q 
h2  k 


Then  is  a  consistent  approximation  to  the  operator  L  described  in  (1.1), 
We  assume  h  is  so  small  that  ak  >  0,  yk  >  0  for  all  k. 

We  are  concerned  with  the  system  of  linear  equations 


(3.4) 


[LhU]k  =  fk,  k  =  1,2 . 2N+1  , 

U0  =  U2N+2  =  0 


vNv\w.v.v.v.v:v:\v 


We  shall  solve  this  system  with  a  particular  multi-grid  iterative  scheme. 
Consider  a  course  grid  on  which  we  have  a  mesh  spacing  of  2h  and  the 
course  grid  points  are  given  by 

Xk(2h)  =  2hk ,  k  =  0,1 ,2 . N+l  . 

We  have  a  space  $h  of  grid  functions  { ( h ) ;  k  =  0,1 , . . . ,2N+2}  defined 
on  the  fine  grid,  and  we  have  a  space  S2h  of  grid  functions 
{ Uk ( 2h ) ;  k  =  0,1,. ...N+l)  defined  on  the  coarse  grid.  Our  first  step  is 
to  construct  the  interpolation  operator  I^h  which  maps  S2h  into  S^, 
i  .e. 

I^1  •  S  -*■  S 
A2h’  *2h  n  * 

Following  the  experience  of  Dendy  [4],  [5]  we  choose  the  following  mapping 
(3.5a)  [l!)hU(2h)]2k  =  Uk^h)  (common  points) 

and 

(3.5b)  [l5hU(2h)J2k.,  •  «,.,(»)  *V2k.,Uk  (2h)]  (new  points). 

This  choice  of  "operator"  interpolation  may  be  described  in  the  following  way 
If  the  physical  point  x.(h)  of  the  fine  grid  is  also  a  physical  point  of 

<J 

the  coarse  grid,  i.e.  if  j  is  even,  say  j  =  2k,  we  set  U ^ ( h )  =  Uk(2h) 
to  be  the  same  value  as  the  coarse  grid  function  assumed  at  that  point; 
that  is,  (3.5a)  holds.  If  the  physical  point  x.(h)  of  the  fine  grid  is 

J 

not  a  point  of  the  coarse  grid,  i.e.  j  is  odd,  say  j  =  2k  -  1 ,  we  re¬ 
quire  that 


(3.5c) 


aji5hu(2h)]),l , 


=  0  . 


We  formalize  this  remark  as 


Lemma  3.1 :  Let  this  interpolation  operator  be  defined  by  (3.5a), 

(3.5b).  Then  a  function  U(h)  e  Sh  is  in  the  range  of  l£h  if  and  only  if 

[Lhu(h)]2k_i  =0.  k  =  1 »2, . . . ,N+1  .  ■ 

2h 

We  now  turn  to  the  construction  of  a  projection  operator  1^  :  Sh  ■* 

We  define 


(3.6)  Cl‘nU(h)]k  -  1/2  U2k.,(h)  -U2k(h)  ♦  U2w(h)  . 


Remark;  if  b(x)  =  0  and  the  operator  L  is  self-adjoint  then 


and  we  see  that 


(3.7) 


ak  *  Vl 


I2h  =  c[Ih  ]T 
lh  CLi2hJ  ’ 


As  we  have  said  in  section  2,  the  relationship  (3.7)  between  1^  and  l"h 
is  the  "variational  choice"  and  is  frequently  recommended  for  self-adjoint 
problems  -  see  [13],  [7],  [4], 

For  the  purposes  of  this  exposition  we  take  the  optimal  choice  of 
"coarse  grid  correction",  i.e. 

i  P  i2hi  rh 

L2h  L2h  ‘  h  Lh x 2h  * 


^  L%  -  "V  m  '•  J,  k  »  "  »  *  ►  *  •  “  ■> 


Remark:  A  direct  computation  shows  that 


(3-8a)  [L2hU(2h)]k  =  -otk  ^Uk-1  +  Sk  \  'Yk  \+l 

where 


(3.8b) 


2ka2k-l 

B2k-1 


(3.8c) 


R(2)=J_fl  a2kY2k-1  Y2ka2k+1 
k  2  2k  B2k_1  ^2k+l 


( 3 . 8d ) 


Y2kY2k+l 

&2k+l 


Hence,  L2h  is  again  a  diagonally  dominant  three  term  operator  of  the 
form  (3.3a). 

Now  we  need  only  describe  the  smoothing  operator  which  we  do  in  the 


next  section. 


4.  Jacobi  Iterative  Schemes 


A  direct  iterative  method  for  the  solution  of  (3.4)  is  described  by 
a  splitting  of  the  operator  .  We  set 

Lh  =  M  -  N  . 

Then,  given  a  first  guess  we  define  successive  iterates  by  the  formula 

(4.1 )  MUJ+1  =  NUj  +  f  . 

The  convergence  of  this  scheme  is  determined  by  the  eigenvalue  problem 

(4.2)  XMU  =  NU  . 

As  is  well  known,  the  scheme  (4.1)  is  convergent  iff 

max|A|  <  1  . 

It  is  easy  to  verify  that:  if  <A,4>>  are  an  eigenpair  of  (4.2)  then 

(4.3)  Lh<{>  =  (l-A)M<t>  . 

In  this  section  we  are  concerned  with  a  particularly  simple  class  of 
such  iterative  methods,  the  Jacobi  methods.  We  set 

(4.4a)  M  *  (l+a)B  ,  a  >_  0  , 

where 

(4.4b)  B  =  diag{9>y)  . 


In  this  case  we  may  rewrite  (4.1)  as 


Uj+1  =  u>j  +  T7a  B_1(f-Lhuj)  • 

When  a  =  0  we  call  this  scheme  the  Jacobi  method.  When  a  >  0  we  call 
this  scheme  a  damped  Jacobi  method.  In  these  cases  we  are  able  to  give 
a  relatively  complete  discussion  of  the  eigenvalue  problem  (4.2).  We  have 
the  following  facts. 

i)  The  method  is  convergent  for  a  _>  0  . 
ii)  Let  a  =  0  .  Let  < y >  be  an  eigenpair,  i.e. 

yfty  =  ( B—  L^ )  cp  . 

Let  $  be  defined  by 

(4.6)  $k  =  (‘1)k+\  * 

Then  <-y,$>  is  also  an  eigenpair.  The  eigenvalues  y  are  real 
and  distinct,  furthermore:  as  h  ■*  0  the  {y >  fill  out  the  interval 
[-1,1].  For  completeness  we  repeat  the  basic  relationship  between 
<p  and  $  .  Namely, 

(4.7a)  if  k  is  odd:  <J>k  =  $k 

(4.7b)  if  k  is  even:  4>k  =  -$k  . 

Since  dim  =  2N+1  there  is  a  single  eigenvector  <j>  associated 
with  the  eigenvalue  y  =  0  .  This  eigenvector  satisfies 


For  a  >  0  there  are  corresponding  eigenpairs  <A,<)>>,  <A,4>>  where 
$  and  <J>  are  the  eigenfunctions  described  in  (2)  and 


(4.9) 


,  _  u+a  t  _  a^i 
A  "  1+a  ’  A  "  1+a 


the  eigenpair  <0,<j>>  corresponds  to  an  eigenpair  (A,<j>>  with 


We  now  turn  to  the  determination  of  the  three  important  subspaces 
Range  l!^»  Nullspaae  Nullspaae  I  j^L^  . 

l_  o  L. 

where  the  operators  1^,  1^  are  given  by  (3.5a),  (3.5b)  and  (3.6). 

Lemma  4.1:  Let  <y,<{)>,  <-y,$>  be  the  two  eigenpairs  described  in 
(ii)  above  with  y  f  0.  Let 

(4.10)  $  =  [(l+y)4>  -  (l-y)$] 

then  i  e Range  .  Further,  since  the  vectors  $  corresponding  to 
different  pairs  <y,4>>,  <-y,£>  are  linearly  independent  this  construc¬ 
tion  provides  N  linearly  independent  elements  of  the  Range  . 

Proof;  Using  (4.3)  we  have 

Lh$  =  [(l+y)(l-u)B4>  -  (l-y)(l+y)B$] 

=  ( 1  +y )  ( 1  -u  )B[4>-$]  . 

Thus,  the  lemma  follows  from  Lemma  3.1  which  characterizes  Range  (i!^) 
and  from  (4.7a),  (4.7b).  | 


Lemma  4.2: 


i)  For  all  choices  of  a  >  0,  the  vector  B<j>  (associated  with  p=0) 
is  an  element  of  Nullspace  1^  and  hence  (4.3)  implies  that  <j> 
is  an  element  of  Nullspace  I^jnL^. 

ii)  Let  <p,<j>>,  <-p,$>  with  p  f  0  be  the  two  eigenpairs  described 
above.  Let 

(4.11a)  4*  =  B [ ( 1  -pi ) 4>  +  ( 1  +y ) <t> ]  . 

Then  ip  is  an  element  of  Nullspace  of  1^  and  (4.3)  implies  that 

(4.11b)  (<J>+$)  e  Nullspace  of  ij^L^  . 

Further,  since  the  vectors  V  associated  with  different  pairs 

<p,<J>>,  '-p,$>  are  linearly  independent,  we  have  N  linearly 

independent  elements  of  this  nullspace  and  N  linearly  independent 

2h  ~ 

vectors  of  Nullspace  1^  L^.  The  vectors  B<{>,  d>  provide  one  more 
independent  vector  of  each  of  these  subspaces. 

Proof:  The  result  follows  from  a  direct  computation  using  (3.6)  and  the 
defining  eigenvalue  problem.  | 

Corollary:  There  is  a  unique  decomposition  of  into  Range  I^h  and 
Nullspace  (ijj^L^.) 

Proof:  Since  Lh  is  non-singular  we  have  shown  that 

dim  Range  N  , 

dim  Nullspace  (ijj^L^)  =  dim  Nullspace  (ij^)  ^  N  +  1  , 
dim  =  2N  +  1  . 

Th"S  the  corollary  follows  from  the  observation  that  Lemma  3.1  implies 
that  these  two  subspaces  have  only  the  zero  vector  in  common.  ■ 


5.  Some  Estimates 


Let  a  ^  0  .  We  take  as  our  smoother  m  applications  of  the  corres¬ 
ponding  Jacobi  iteration.  That  is,  given  U°  =  U°(h)  we  obtain  0 
(as  in  2.6b)  from  the  formula 

(5.1a)  MUj+1  =  NUJ  +  f  ,  j  =  0,1 .... ,m-l  , 

(5.1b)  0  =  1/". 

First  let  us  consider  the  special  vector  $  with  its  associated  eigenvalue 


X 


Suppose 


0 

e 


then 

e  =  CXm$ 
lhe  =  CXmB$ 

and,  using  Lemma  4.2  we  see  that 

r  2  h  -  ~  ^ 

!h  V  ■ 0  • 

Hence,  in  this  case,  for  any  norm 


ft 


Me  now  consider  two  norms  defined  on  .  Let  w,  v  e  '  be  of  the 


(5.3a) 


_  it  A 

w  =  l  A  .  +  A4>  +  l  A  . 
j=1  J  J  j=l  J  J 


(5.3b) 


V  =  l  C.A.  +  +  l  C.$.  . 

j=1  J  J  j=l  J  J 


Define 


N  ~~  N  ^ 

<w,v  >0  *  l  A  C  +  AC  +  l  A  C 
U  j=l  J  J  j=l  J  J 


<w,v>  =  l  A.C.(l-U-)  +  AC  +  l  A.C.(l+y.) 

1  j=l  J  J  J  j=l  J  J  J 


w||0  =<w,w>0 


=  <w,w>,  s  l  A^(l-y . )  +  A2  +  l  A2(l+y.) 
i  j=l  J  J  j=1  J  J 


Lemma  5.1:  Suppose 


(5.4) 


on  a 

€u  =  u  -  U  =  C4>  +  d<f> 


and  let  £  be  defined  by  the  two-grid  iteration  scheme.  Then 


(5.5a) 


II  e1 1|  2 


max - n — 

C.d  ||  e°||g 


0  _ 


=  1/2 


2m  2  2m  2 

O  +  ^  (1+u) 


‘  1+a' 


(5.5b) 


max 

c.d 


1 II2 

e  Hi 

“U7T 

e  lh 


=  1/2 


02md-.)  +  (f^O+n) 


Proof;  From  (5.4)  we  see  that 


e  =  cXm<{)  +  d^$  . 


Following  the  theory  of  section  2,  we  write  (see  2.16) 
(5.6) 


?  =  l2hw](2h)  +  £](h)  . 


where 


We  claim  that 
(5.7a) 

(5.7b) 


'hV'h>  ■  ° 


l5h»’(2M  ■  [(HpH-(l-u)?] 


U4]  . 


To  verify  this  we  need  merely  verify  that  the  sum  of  the  right-hand-sides 
of  (5.7a)  and  (5.7b)  is  ?  ,  and  use  (4.10)  and  (4.11b). 


Having  verified  (5.6),  (5.7a),  (5.7b)  we  proceed  as  follows 

(5.8a)  ||  £]||q  =  l/2[cAm(l-y)  +dAm(l+y)]2 

(5.8b)  ||  e°||  jj  =  c2  +  d2 

(5.9a)  ||  e1||^  =  l/2[cXm(l-y)  +dXm(l+y)]2 

(5.9b)  ||  €°||  2  =  C2(l-y)  +  d2(l+y) 

Thus,  a  simple  argument  shows  that 

(5.10a)  sup  1^0  .  1/2[x2m(l-y)2  n-X^d+y)2]  , 

II  €  Ho 

m  1 1|  2 

(5.10b)  sup  \  l  =  l/2[X2m(l-y)  +X2m(l+y)]  . 

II  «°llf 

Using  the  basic  formulae  (4.9)  we  obtain 


sup 


sup 


Iloilo 

M  -1  II  f 

iT?iiT 


=  1/2 


=  1/2 


OaB(1_u)2  +(^)2md+y)2 


W 


(ir|)2m{1‘y)  +(4ir)2m(i+y) 


T+a' 


Thus,  the  lemma  is  proven.  I 


Theorem  5.1:  In  the  general  case  we  have 


(5.11a) 


<  2UP  •  1/2  <ij|)2l"n-u)2  +(yjj)2m(l*u)Z 


(5.11b) 


lFT®ih  1  -i£i  ,/2P  (Uu)  t(^,a"(,tu) 


Proof:  The  Theorem  follows  immediately  from  the  previous  lemma.  I 

We  observe  that  (5.11a)  is  precisely  the  formula  obtained  by  Hackbusch 
[12,  (2.21)]  in  the  special  case  p(x)  =  1,  b(x)  =  q(x)  =  0  and  a  =  1  . 
To  make  a  complete  identification  we  merely  set 


(5.12) 


O  =  Hi,  (Uo)  .Hi  . 


However,  while  (5.11a),  (5.11b)  describe  the  worst  case  decay  in  one  multi¬ 
grid  iteration  in  a  2-grid  scheme,  it  does  not  give  the  estimate  of  real 
interest.  From  the  discussion  in  the  proof  of  Lemma  5.1  -  and  (5.7)  in 

particular  -  we  realize  that,  even  though  the  constants  c  and  d  of  (5.4) 

0  k 

may  be  arbitrary  for  e  ,  that  is  not  true  for  e  ,  k  >  1  .  We  have 
(from  (5.7b)) 


_  _  cA  (1-m)  +  dA^'n  +u) 

a  "  2 


Therefore,  following  the  argument  of  Lemma  5.1, 


e{2)  =  f[Am(l-u)  +  Am(l+y)]  [<£+<{>] 


Hence 


IU(,)|l^  *  l|e(”ll?  -  2c2 


l|e(2)||2  =  I|e(2)||2  =^[X"(1-U)  ♦fOn.)]2 


and  for  j  =  0,1  we  have 


-(2)  II  2 


4  =  jtAl-y)  +An,(l+y)]2  . 


(5.14) 


Thus  we  have  proven 


Theorem  5.2;  In  the’  general  case,  for  j  =  0,1  and  all  k  >  1  we  have 


.(k+1 )  i,  2 


(5.15) 


|[xm(l-y)  + Am 


(i+y)]2}  . 


Remark:  The  distinction  between  Hackbusch's  result  (5.11a)  and  (5.15)  is 
non-trivial  for  large  m.  We  have,  as  m  -*•  ® 


I  C 


'o  _  (1+a)  1 


6.  Symmetrization 

Consider  the  difference  equation  (3.4)  described  by  (3.3a),  (3.3b) 


Let 


(6.1) 


=  dkVk 


where  the  coefficients  d^  are  computed  recursively  by 


(6.2a) 


=  1 


(6.2b) 


k=0,l .... 


then  (3.4)  becomes 


"akdk-l Vk-1 


+  6kdkVk  '  Ykdk+lVk+l  =  fk 


or 


(6.3) 


-a 


k 


+  ¥k 


k+1 


k+1 


that  is 


26 


where 

(6.8b)  D_1f  =  (?!  ,?2 . ^2N+-|)T  * 

Since  B  is  also  the  diagonal  of  Lh>  (6.8)  is  precisely  the  same 
Jacobi  iterative  scheme  for  modified  symmetric  equations. 

Now  let  us  study  the  effect  of  this  transformation  on  the  two-grid 
iterative  scheme.  We  compute 

Dk^I2hU^k  * 

Imagine  U ( 2h ) k  =  d2kV ( 2h ) k  is  given  in  S2h  .  Then,  from  (3.5a) 

(6.9a)  d2k[I2hU(2h)]2k  =  V(2h)k  * 

(6.9b)  d2k-2^2hU^2h^2k-2  =  V^2h^k-1  * 

Thus  (3.5b)  yields 

d2k-1^2hU(2h^2k-l  =  6^~j‘^ct2k-ld2k-ld2k-2Vk-l)  +Y2k-ld2k-ld2kVk 


=  6^7r[ot2k-lV(2h)k-l  +Y2k-lV(2h)k^  * 


Thus,  with  this  change  of  variables  the  mapping  1^  of  our  original 

?h 


unsymmetric  problem  becomes  I2h»  the  appropriate  mapping  associated 
with  the  new  symmetric  problem. 


Finally,  let  us  consider 


d^Ci2hhu(M]k  • 

A  straight  forward  calculation  verifies  that 


d2k^h  U-*k  =  2^8 


2h 


2k-l 


V2k-1^  +v2k(h)  +  8 


^2h 

2k+l 


V2k+1-* 


Thus,  following  the  remarks  of  section  3  [see  (3.7)]  we  see  that 


l-H-Vhi 
Cd2kIh  ] 


9 


the  appropriate  projection  operator. 

For  our  purposes,  the  major  significance  of  these  calculations  is  that 
the  "1"  norm  introduced  in  section  5  is  the  "operator  norm"  for  the 
symmetric  problem.  Hence,  we  have  a  norm  which  is  well-defined  on  all 
spaces  Sh  . 


Multi-Grid  and  Experimental  Results 


The  results  of  the  previous  sections,  and  Theorem  5.1  in  particular, 
provide  exact  estimates  of  the  decay  of  the  error  (in  two  norms)  in  one 
iteration  of  a  2-grid  scheme  -  in  the  worst  case. 

Since  l_2^  is  again  a  three  term  (diagonally  dominant)  operator  of 
the  form  (3.3a)  -  and  given  specifically  by  (3.8)  -  we  may  apply  our 
multi-grid  approach  inductively  as  follows:  Assume  that  the  n-grid  multi¬ 
grid  scheme  based  on  "smoothing"  with  m  applications  of  the  damped  Jacobi 
iteration  with  parameter  a  is  defined.  Suppose 


(7.1a) 


h  =  —  ,  n  >  2  , 
2n 


where  H  is  of  the  form 


(7.1b) 


H  =  pTT  »  P  1  1  ' 


We  wish  to  solve  (3.4)  on  the  h-grid.  The  iterative  scheme  is  given  by  the 
following  inductive  description. 

(1.)  On  the  h-grid  (h  =2~nH): 

(a)  Let  be  chosen.- 

(b)  Apply  the  damped  Jacobi  (with  parameter  a)  iteration 
m  times  to  obtain  U. 

(c)  Form  r(h)  =  f  -  LhU. 

(2.)  Transfer  Information: 

(a)  Set  r(2h)  =  ljjhr(h) 

(3.)  On  the  2h  grid: 

(a)  Consider  the  problem 


W 


(b)  if  2h  =  H,  solve  exactly. 

(c)  if  2h  <  H,  set  U°(2h)  =  0  and  apply  the 
n-grid  iterative  scheme  (based  on  smoothing 
with  m  applications  of  the  damped  Jacobi 
iteration  with  parameter  a).  Let  U^(2h) 
be  the  result  of  this  step. 


(4.)  Transfer  Information: 

(a)  U1  =  0  +  l!jhu1(2h)  . 

(b)  U1  -  U°  . 

Return  to  1(b). 

In  the  multi-grid  jargon  this  is  the  so-called  slash  or  sawtooth  cycle 
which  we  indicate  schematically  as: 


/  h 
2h 

ST  4h 


h  o 

\  m  "smoothings" 

2h 

m  "smoothings" 
4h 

°s 

2H  sa 

pi - 1 

H  solve 


Note:  There  are  no  smoothing  steps  during  the  transfers  from  coarse  to  fine 
grids. 

McCormick  [14],  [15],  calls  such  a  multigrid  cycle  a  Mv h  cycle.  When 
the  smoothing  occurs  only  on  the  way  "up"  the  cycle  and  the  errors  are 
merely  restricted  on  the  fine  to  coarse  transfers,  he  calls  the  cycle  a  M /. 


cycle.  For  the  symmetric  case,  using  Richardson's  iteration,  see  [14], 
he  shows  that 


('•2)  l|H/hll,  -  l|Mxhl|,  . 

In  discussing  the  symmetric  cycle  he  obtains  the  following  estimate. 
Let 


(7.3a) 


0 

£ 


0 

n 


+  I 


h 
2ha) 


Suppose  a,  0  <  a  <  1  satisfies 

(7.3b)  !|Ge°llf  ^  c.||n°||f  +  ||  I^hu>°||f  .  Ve°, 


then 


(7.4a)  11%!!,  <aH, 

that  is: 

(7.4b)  lie'll,  »  l|u-U'||,  <  ot**  ||  e°||,  . 

Since  our  Jacobi  iteration  is  not  all  that  different  than  Richardson's 
iteration  it  is  not  surprising  that  a  similar  result  holds  in  our  case. 
Indeed,  if  one  applies  his  argument  to  our  multigrid  cycle,  i.e.  the  M^h 
cycle,  one  gets  the  following  result. 


Lemma  7.1:  Consider  the  symmetric  case  and  suppose 

(7-5)  II  G|l7  <  1  - 

Let 

(7.6a)  Ge^  =  e  =  n  +  . 

2h 


Suppose  a,  0  <  a  <  1  satisfies 


(7.6b)  II  n II f  ♦  aili^lf  <  S||  c°||f. 

Then 

(7.7)  ||Mxh||i  <  S*5  • 

Proof;  The  proof  proceeds  by  induction.  Since  we  use  an  exact  solver 
on  the  coarsest  grid, 

(7.8a)  II  M\h  ||  i  *  0  <  S’5  . 

Assume  that 

(7.8b)  ilMx2hn1  ^  a*5 , 

that  is: 

(7.8c)  ||  MN2hw(Zh)-u(2h)||1  <  a*||w(2h)  ||  1  . 

Then 

(7.9)  II  e1  ||f  =  ||  n  ||f  +  II  I2h(M\2hw“u,)Hf  - 

Because  this  is  a  symmetric  problem  we  know  that 
(7.10a)  ||  v  || f  -  <v,Lv>  , 


and  that 


32 


(7.10b)  II  1 2h» II f  t-hI2h',>  =  <v-  ^Ih''LhI5h''> 

■  Y<».  4h',)  ■  rIMI*  • 

Therefore, 

lie’ll?  -  II  n  ||  f  ♦  ill  M\2h“"“lll  • 

By  the  inductive  hypotheses  (7.8b)  we  have 

(7.11)  He’ll*  <  Hull?  ♦  v S II o || *  *  II n || ?  +  S||ljh«i||* . 

Note:  In  (7.10b)  and  in  this  calculation  the  symbols  ||u||^  and  il  ^ 2hu ^  1 
refer  to  the  designated  norms  on  the  spaces  S2h  and  . 

By  the  basic  inequality  (7.6b)  we  have 

lie’ll,2  <  || n || 2  +  “II l|h “ II ?  i  S||  e°||? 

which  proves  the  Lemma. 

Since  the  proof  of  this  lemma  is  immediate  once  one  understands  the 
proof  of  McCormick's  lemma  2.2  of  [15]  one  would  expect  that 

(7.12)  a  -  a,  . 


Indeed,  this  is  the  case.  Direct  but  messy  calculations  based  on  the  results 
of  section  5  yield 


(7.13)  a  =  a  = 


sup 

-1<U<1 


,  2m  ,  ^  2m  2m  2m 

0-p)-O 

(uf] 


1  -  n*M) 


Moreover,  for  all 

choices  of 

a  and 

m,  the 

su premum  is 

attained  at  p  =  1  . 

i* 

The  corresponding  values  of  a?  are 

displayed 

in  the  following  table: 

Bounds  on 

the  Convergence 

Rate 

m\a 

.333 

.5 

.667 

.75 

1  1 

.333 

1 

.633 

.577 

.561 

.561 

.577 

.614 

2 

.435 

.408 

.417 

.424 

.447 

.475 

3 

.336 

.335 

.349 

.357 

.378 

.403 

4 

.283 

.293 

.307 

.314 

.333 

.357 

In  view 

of  the 

results  of 

section 

6  which 

demonstrate 

the  complete  equiva- 

lence  of  our  problem  to  a  related  symmetric  problem,  these  upper  bounds  apply 


in  our  case. 

\c  /s  i,c 

However  the  estimate  =  a  ^  is  only  an  upper  bound  for  the  rate  of 
convergence  of  the  multigrid  iterative  scheme.  In  order  to  complete  our  in¬ 
vestigation  we  have  undertaken  an  experimental  project. 

A  computer  program  was  written  with  the  following  capabilities:  The  user 
suppl ies 

p(x),  b(x) ,  q(x) ,  f(x),  m,  a,  n,  and 

M  =  (-^)  -  1  =  (number  of  points  on  the  finest  grid), 

where  p(x),  b(x),  q(x),  f(x)  are  the  coefficients  of  the  problems  (1.1), 
(1.2)  and 

m  =  number  of  applications  of  the  damped  Jacobi  iteration, 
a  =  parameter  of  the  damped  Jacobi  iteration 
n  =  #  of  grid  levels. 

The  user  also  supplies  an  initial  guess  U°  and  a  tolerance  E. 


The  program  then  executes  multi-grid  iterations  until  the  norm 
[see  (7.14a)]  of  the  residual  is  below  the  indicated  tolerance  E  .  The 
program  is  run  in  an  interactive  fashion  which  allows  the  user  to  change 
the  parameters  M,  m,  a  and  n  . 

The  experiments  reported  here  were  run  on  the  VAX  780  in  both  single 
and  double  precision  arithmetic  (approximately  sixteen  decimal  digits  of 
accuracy).  The  single  precision  results  were  qualitatively  similar  to  the 
double  precision  results,  however,  for  increased  accuracy,  the  double  pre¬ 
cision  results  are  reported  here. 

For  our  present  purposes  the  basic  program  was  modified  to  enable  us 
to  estimate  the  "rates  of  convergence"  of  the  multi-grid  iteration.  For 
each  test  problem  we  used  a  known  solution  u(x)  of  the  boundary  value 
problem  (1.1),  (1.2).  Then  we  computed  the  exact  solution  u(x,h)  of 
the  algebraic  system  (2.4).  Then  using  two  norms 

(7.14a)  ||u||  =  hi  |uj| 

(7.14b)  ||  u||  ]  =  <  u»Lhu  >** 

we  computed  the  norms  of  the  error  (u-u1)  at  each  iteration.  The  rate  of 
convergence  was  measured  by  computing 


at  each  iteration  i  =  1,2,3.... 

To  check  that  the  program  was  working  correctly  a  number  of  measures 
were  taken.  The  most  simplistic  was  to  carry  out  some  of  the  iterations 
by  hand  and  to  compare  the  hand  computed  calculations  to  the  iterates 


generated  by  the  machine.  In  addition,  since  the  discretization  error 
2 

is  0(h  ),  it  is  not  unreasonable  to  expect  that  halving  the  step  size 


should  reduce  the  final  error  in  u  by  a  factor  of  four.  This  property 

2£+l  h 

was  checked  and  found  to  be  true.  One  of  the  requirements  for  I  0 

2  n 

is  that 

<7->6>  |lai^,“  - 0 

2  h  2  h  J2k-1 


( by  1 emma  2.1). 

After  each  coarse  to  fine  grid  transfer,  formula  (7.16)  was  computed  and 
checked.  Finally,  from  (5.7b)  one  sees  that  the  error,  e2k(h),  on  tfle 
even  points  of  the  coarsest  grid  should  be  zero.  This  requirement  was 
also  verified  after  each  iteration. 

The  test  problems  are  best  described  by  giving  the  choices  of  p(x), 
b(x),  q(x)  and  u(x),  the  true  solution  of  the  differential  equation 
(1.1),  (1.2)  (which  determines  f(x)). 

As  a  basic  case  we  took 


(7.17a)  p(x)  =  1,  q(x)  =  b(x)  =  0  and  u(x)  =  0. 

This  test  was  merely  to  be  sure  the  program  worked  on  this  simple  case. 
In  addition  there  were  six  other  problems  based  on  two  additional  sets 
of  coefficients  p(x),  b(x),  q(x)  and  three  "solutions"  u(x).  These 
are 

1  2 
(7.17b)  p ( x )  =  1  +  2  sin  4irx,  b(x)  =  1  +  x,  q(x)  =  (sin  5ttx) 


(7.17c)  p(x)  =  ex,  b(x)  =  1  +  x2,  q(x)  =  (l-x)eX/>2 


The  "solutions"  were 


(7.18a) 

u^x)  =  x(e-ex)  , 

(7.18b) 

u2(x)  =  x5/2(l-x)  , 

(7.18c) 

u3(x)  =  sin(14Trx)  . 

For  each  problem,  test  runs  were  made  with  a  variety  of  initial  guesses 
After  all,  the  point  was  to  obtain  the  worst  rate  of  convergence.  Each 
initial  guess  consisted  of  a  smooth  component 

S  kTT 

u.  ■  20  sin  W7T-  where  M  is  the  number  of  points  on 

the  finest  grid 

and  a  rough  component.  The  rough  component  was  chosen  in  various  ways 
in  order  to  have  different  compositions  on  the  coarser  grids.  The  rough 
components  of  the  initial  guesses  are  best  described  schematically,  by 
setting 

\  *  uk  4  405k 

where  |  6^1  =1,  and  the  sign  of  <$k  follows  the  following  patterns: 

Initial  Guess  Pattern  for  6 
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Runs  were  made  with  a,  the  damped  Jacobi  parameter,  equal  to 
.333,  .5,  .667,  .75,  1.0,  1.333,  while  m,  the  number  of  smoothing 
iterations,  ran  from  one  to  four  and  the  number  of  grid  layers  varied  from 
two  to  five.  For  each  test  problem,  the  program  stopped  when  the  discrete 
norm  of  the  residual  vector  was  less  than  .00005.  The  most  recently 
computed  rate  of  convergence. 

II  e final  H 
II  efinal-l  11^ 

was  computed  and  recorded  in  Tables  I I I- VI . 

The  theoretical  rate  for  a  two  grid  iteration  scheme  was  computed  from 
Theorem  5.1  and  Theorem  5.2  by  solving  for  the  maximum  of 


and 


F(V) 


1 

2 


u+a] 

1+aJ 


2m 

(l-y) 


+ 


+ 


-1  <  y  <  1 


-1  <  y  <  1 


using  Newton's  method.  Table  I  exhibits  (max  F(y)P  (a  predicted  rate  of 
convergence)  as  a  function  of  m  and  a  .  The  value  of  y  at  which  the 
maximum  of  F(y)  occurred  can  be  found  in  Table  I' 
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Table  I 


Predicted  Rate  Based  on  F(y) 


m\a 

.333 

.500 

.667 

.750 

1 .000 

1.333 

1 

.500 

.333 

.400 

.429 

.500 

.571 

2 

.260 

.248 

.261 

.268 

.289 

.331 

3 

.200 

.206 

.217 

.223 

.238 

.258 

4 

.171 

.180 

.190 

.195 

.208 

.225 

Table  I* 

m\a 

.333 

Damped  Jacobi 

.500  .667 

Parameter- 

.750 

y 

1.000 

1.333 

1 

1 

0 

0 

0 

0 

0 

2 

.883 

.707 

.666 

.650 

.577 

.370 

3 

.833 

.786 

.762 

.750 

.714 

.661 

4 

.857 

.833 

.814 

.811 

.777 

.741 

Table  II  exhibits  (max  F-^y))1*  (a  better  rate  of  convergence)  as  a  function 
of  m  and  a.  The  value  of  y  at  which  the  maximum  of  F.j(y)  occurred 
can  be  found  in  Table  II*. 

Table  II 


Predicted  Rate  Based  on  F^(y) 


333 

.500 

.667 

.750 

1.000 

1.333 

500 

.333 

.400 

.429 

.500 

.572 

250 

.111 

.160 

.184 

.250 

.326 

125 

.078 

.088 

.093 

.125 

.187 

,068 

.062 

.068 

.072 

.083 

.109 

Table 

III 

Worst 

case,  2-grids 

m\a 

.333 

.5 

.667 

1 

.  499 ( 3 ) 

,333*b* 

.400^ 

2 

.250(c) 

.160(b) 

3 

.124(c) 

.075<»> 

. 087  ^ d  ^ 

4 

.063(a) 

■062*f* 

.068(g^ 

Table 

JV 

Worst 

case,  3-grids 

m\a 

.333 

.5 

.667 

1 

.499(c) 

.333{b) 

.400<b> 

2 

.250(c) 

,165(9> 

.192(g) 

3 

.124(c) 

.099^ 

.1 16(k) 

4 

.087(j) 

.079^ 

.087^ 

Table 

Worst 

case,  4-grids 

[( 


Worst  case,  5-grids 


m\a 

.333 

.5 

.667 

.75 

1.0 

1.333 

1 

.499(0 

.333(b) 

.400(b) 

.429(b) 

.500(b) 

.571  ^ 

2 

.250^ 

.209(o) 

.221 ^ 

.222^ 

.268(k) 

.337(g^ 

3 

,124^c) 

.134(o) 

.148(o) 

.148^ 

.175(k) 

2i9(k) 

4 

.098^ 

.098(p) 

.104^ 

.lll(k) 

.131{k) 

.160(k) 

The  letters  in  the  above  tables  correspond  to  the  choices  of  coefficients, 
"solutions",  and  patterns  for  rough  components  in  the  initial  guess  [see 
(7.17),  (7.18),  (7.19)]  displayed  in  Table  VII. 

Concluding  Remarks 

As  can  be  seen  from  the  computational  results,  no  particular  choice  of 
problem  or  initial  guess  always  resulted  in  giving  the  worst  case.  Moreover, 
it  appears  that  a  =  a  is  aji  upper  bound  on  the  rate  of  convergence  of  the 
multigrid  scheme  but  does  not  yield  the  exact  rate  of  convergence.  Notice 
that  there  seems  to  be  no  degradation  for  m  =  1.  However,  as  m  increases 
we  find  some  degradation  in  the  rate  of  convergence.  But,  it  appears  to  be 

i _ 4. i _  *5  _  C  *5 


Table  VII 


Worst  Case  Problems 


Problem 

Coefficients 

"Solution" 

Pattern  for 
in  initial  guess 

a 

7.17b 

x(e-ex) 

B 

b 

7.17c 

sin(14irx) 

B 

c 

7.17b 

5/2, ,  » 

x  (1-x) 

C 

d 

7.17c 

sin(14Trx) 

C 

e 

7.17b 

x{e-ex) 

C 

f 

7.17c 

x(e-ex) 

D 

9 

7.17a 

0 

D 

h 

7.17c 

x(e-ex) 

D 

i 

7.17c 

sin(14TTx) 

E 

j 

7.17a 

0 

E 

k 

7.17b 

x(e-ex) 

E 

l 

7.17c 

x5/2(l-x) 

A 

m 

7.17a 

0 

E 

n 

7.17c 

x5/2(l-x) 

C 

0 

7.17b 

x(e-ex) 

A 

P 

7.17c 

x(e-ex) 

A 

q 

7.17b 

sin(14irx) 

A 

'  »**  h'  i *»  A-'  v' 
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